An Adjoint Approach for Stabilizing the Parareal Method

The parareal algorithm seeks to extract parallelism in the time-integration direction of time-dependent differential equations. While it has been applied with success to a wide range of problems, it suffers from some stability issues when applied to non-dissipative problems. We express the method through an iteration matrix and show that the problematic behavior is related to the non-normal structure of the iteration matrix. To enforce monotone convergence we propose an adjoint parareal algorithm, accelerated by the Conjugate Gradient Method. Numerical experiments confirm the stability and suggest directions for further improving the performance. To cite this article: F. Chen, J.S. Hesthaven, Y. Maday, A. Nielsen, C. R. Acad. Sci. Paris, Ser. I ??? (2015). Résumé L’algorithme parareel en temps vise à proposer la parallélisation dans la direction temporelle de l’intégration approchée déquations différentielles instationnaires. Cet algorithme a été appliqué avec succès sur une large gamme de problèmes mais souffre néanmoins de certains problèmes de stabilité lorsqu’il est appliqué à des problèmes non dissipatifs. Nous exprimons l’algorithme comme une méthode d’itération matricielle et nous montrons que le comportement problématique est liée à la structure non-normale de la matrice d’itération. Pour retrouver une convergence monotone nous proposons une version adjointe de l’ algorithme pararéel, accélérée par la méthode du gradient conjugué. Des expériences numériques confirment la stabilité et proposent des pistes pour améliorer les performances. Méthode adjointe pour la stabilisation de l’algorithme pararéel. Pour citer cet article : F. Chen, J.S. Hesthaven, Y. Maday, A. Nielsen, C. R. Acad. Sci. Paris, Ser. I ? ? ? (2015). Email addresses: [email protected] (Feng Chen), [email protected] (Jan S Hesthaven), [email protected] (Yvon Maday), [email protected] (Allan Nielsen). Preprint submitted to the Académie des sciences August 8, 2015 1. Parareal in time algorithm The increasing number of processors available at large computers presents a substantial challenge for the development of applications that scale well. One potential path to increased parallelism is the development of parallel time-integration methods. Once some standard methods-of-lines approach have been applied, the problem is traditionally viewed as a strongly sequential process. Attempts to extract parallelism have nevertheless been explored : a simple example is parallel Runge-Kutta methods where independent stages allow for the introduction of parallelism [8]. The parareal method, first proposed in [4], is another and more elaborate approach to achieve a similar goal. This algorithm borrows ideas from spatial domain decomposition to construct an iterative approach for solving the temporal problem in a parallel global-in-time approach. To present the method, consider the problem  ∂u ∂t +A (t,u) = 0 u (T0) = u0 t ∈ [T0, T ] (1) where A : R× V → V ′ is a general operator depending on u : Ω×R → V with V being a Hilbert space and V ′ its dual. Assume the existence of a unique solution u (t) to (1) and decompose the time domain of interest into N individual time slices T0 < T1 < · · · < TN−1 < TN = T , where Tn = n∆T . Now, for any n ∈ N, we define the numerical solution operator Fn ∆T that advances the solution from Tn to Tn+1 as Fn ∆T (u (Tn)) ≈ u (Tn+1) (2) This allows, starting from U0 = u0 to define recursively approximations U1, · · · ,UN to u(T1), · · · ,u(TN ) by setting, for n = 0, . . . , N − 1 : Un+1 = Fn ∆T (Un). It is pertinent, for what follows, to realize that this numerical solution corresponds to the forward substitution solution method applied to the matricial system MFŪ = Ū0 where MF , Ū and Ū0 are defined as follows MF =  1 −F0 ∆T . . . . . . . . . −FN−1 ∆T 1  , Ū =  U0 .. .. UN  , Ū0 =  u0 0 .. 0  . (3) If we instead solve the system using a point-iterative approach, i.e., seeking the solution on form Ū = Ū + ( Ū0 −MFŪ ) , we may observe that Ū is known at each iteration, allowing a complete decoupled computation of F1 ∆T , · · · ,F TN ∆T on all intervals. However, for any k, we clearly have to wait k iterations in order to get Uk = Uk hence we have to wait N iterations of the above iterative approach to solve the problem. Consequently, the above approach does not provide any reduction in the time solution. In order to reduce the number of iterations K needed to reach convergence, we need to find an appropriate preconditioner to accelerate the iteration, a typical approach is to utilize an approximation MG ≈ MF , where MG is cheap to apply. We can readily create such an MG by defining another operator G∆T that proposes a coarser approximation Gn ∆T (u (Tn)) ≈ u (Tn+1) (4) by reducing accuracy and choosing a coarse grained model or an entirely different numerical model. Solving the system using standard preconditioned Richardson iterations we may write Ū = Ū + (MG) −1 ( Ū0 −MFŪ ) . (5) This recovers the parareal algorithm on its standard form U n+1 = G Tn ∆TU k+1 n + F Tn ∆TU k n − G Tn ∆TU k n with U 0 n+1 = G Tn ∆TU 0 n and U k 0 = u(T0). (6) 2 This approach has no inherent limits to the amount of parallelism that can be extracted. Important contributions to the analysis of the method can be found in [6,1,3].